Market Location Analysis

4.1.Shapefiles

#required Libraries 
library(sf)
library(sp)
library(sf)
library(raster)
library(dplyr)
library(spData)
library(knitr)
library(googleway)
library(tmap)
library(kableExtra)
library(ggplot2)
library(DT)
library(data.table)

Dataset1: Grid of Norway map consisting of cells where each cell is 1 x 1 km^2 downloaded from statistics bureau of Norway. The dataset Can be downloaded from this link.

ssb1km <- st_read(dsn = "~/Downloads/nhh/spatial data analysis/project/ssb001k" , layer = "ssb1km")

Reading layer `ssb1km' from data source `/home/sherif/Downloads/nhh/spatial data analysis/project/ssb001k' using driver `ESRI Shapefile'
Simple feature collection with 533918 features and 6 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -106000 ymin: 6419000 xmax: 1131000 ymax: 7971000
epsg (SRID):    32633
proj4string:    +proj=utm +zone=33 +datum=WGS84 +units=m +no_defs

The variable SSBID is the ID for each cell in the grid.

Dataset 2: Since the the previous file does not have names of counties or administrative boundaries of norway, and becuase the focus of this project is on Bergen area, the administrative boundaries geo-data of Norway has been downloaded from www.gadm.org.

adm <- read_sf( dsn = "~/Downloads/nhh/spatial data analysis/project/NOR_adm" , layer = "NOR_adm2")

Checking Administrative boundaries data file.

DT::datatable(adm)

Since Bergen area is the one of interest, it has to be separated from the data.

Bergen = adm %>% filter(NAME_2 == "Bergen")

Setting Bergen administrative boundaries data to be the same CRS as the grid data.

Bergen <- st_transform(Bergen, crs="+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs")

Next, the administrative boundaries data of Bergen and the grid data, had to be joined together, so we can have the grid cells of Bergen only.

#Intersect layers
ssb1kmBergen<- st_intersection(ssb1km, Bergen)

Plotting the results.

#plot ssbid
tmap_mode("plot")
tm_shape(ssb1kmBergen) + tm_polygons() + 
tm_credits("Bergen Grid", position=c("right", "bottom"),size = 1.5)

Dataset 3: The total number of population for each cell can be fetched from here.

#\Reading population data 
population = read.csv("Ruter1000m_beflandet_2018.csv" , sep=";", header=TRUE)

Checking population data file.

knitr::kable(head(population,n = 5))

#DT::datatable(population,n = 5)

Variables of interest in this dataset are; SSBID which is the ID for each grid’s cell and “pop_tot” which is the total population in each cell.

To Merge the population data with the grid data, some joining operations had to be done.

#changing the class of cells' ID in both to be charachter 

population$ssbid_1000m= as.character(population$ssbid_1000m)
ssb1kmBergen$SSBID= as.character(ssb1kmBergen$SSBID)

#Joining population data with the shapefile
Bergen_popluation = left_join(ssb1kmBergen, population, by = c("SSBID"= "ssbid_1000m"))

So now, we have a final file called “Bergen_populaiton”, We can visualize its information as follows:

tmap_mode("view")
tm_shape(Bergen_popluation)+ tm_fill(col ="pop_tot",breaks = c(0, 1 ,1000,2000,3000,4000, 5000, 6000, 7000, 8000))

#+  tm_credits("Bergen population for each cell", position=c("right", "bottom"),size = 1.5)

As we see from the map, there are some missing numbers of total population for some grid cells, however, my interpretation that these are unpopulated areas so they can be converted to be zero instead of NA.

Bergen_popluation$pop_tot[is.na(Bergen_popluation$pop_tot)] <- 0

4.2.Bunnpris Geo- coded locations using Google API

After finishing with having geo-data of Bergen area, the next step was to find a good brand candidate as a retailer for this project. “Bunnpris” retail chain has been chosen since it has a reasonable number of branches across Bergen, not so few so we can not have any insights nor too many so it can complicate the mapping.

For retreving the geo-location of Bunnrpis’ branches across Bergen, Google Places API had to be used. For using this service, an API key has to be generated for the use. Instructions for getting the api can be found on this page.

For handling the google service through R interface, googleway library was installed. The following steps were done:

1- Using set_key to make API keys available for all the google_ functions, so there will be no need to specify the key parameter within those functions

# set_key(key = "xxx xxx xxx)

2- Searching for Bunnpris near Bergen area

# bunnpris <- google_places(search_string = "bunnpris near Bergen")

Searching by the string “near bergen”" gives better results than “bergen,hordaland” or“,bergen”.

This function returned a nested list which one of its components is a dataframe contains addresses, geolocations of bunnpris market in Bergen.

Googleway Library also allows for searching with a specific radius around certain centriod.

#bunnpris <- google_places(location = c(60.3913, 5.3221), keyword = "bunnpris", radius = 5000)

Howeve, this function did not give better results than the previous one.

3- As aforementioned, the output of the search string is a multiple nested list. Therefore, a series of operations had to be done in order to get the geocoded locations separately.

#bunnpris.results = bunnpris$results

#bunnpris.geometery= bunnpris.results$geometry
#bunnpris.geocode = bunnpris.geometery$location
#bunnpris.geocode$adress = bunnpris.results$formatted_address
## adding the adress column 

Now, we have a dataframe contains the address and geo-coded location of 18 Bunnpris shops on the area of Bergen. According to Bunnpris website,there are only 16 shops on Bergen area may be they have different definitions of Bergen administrative frontiers.

bunnpris.geocode = read.csv("bunnpris.geocode.csv")
knitr::kable(head(bunnpris.geocode,n = 5))

4- Converting the dataframe into a spatial object.

bunnpris_sf_longlat <- bunnpris.geocode %>% 
   st_as_sf(coords = c('lng', 'lat'), crs = "+proj=longlat +datum=WGS84")



#converting the CRS to the one similar to our base data

bunnpris_sf <- st_transform(bunnpris_sf_longlat, crs="+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs")

Now we can see on the map where are Bunnpris branches.

tmap_mode("view")
tm_shape(Bergen_popluation)+ tm_fill(col ="pop_tot",breaks = c(0, 1 ,1000,2000,3000,4000, 5000, 6000, 7000, 8000)) + tm_shape(bunnpris_sf) + tm_symbols(shape=20,size =.1,col = "green" )  + tm_legend(outside = TRUE, outside.position = "bottom", stack = "horizontal")

5.1.Proximity Analysis

In R there are mainly two packages found to implement Delaunay triangulation, “deldir” and “Dismo”. A comparison between the two packages is done and the output is used for answering the research question.

 coor=as.data.frame(st_coordinates(bunnpris_sf_longlat))
#st_coordinates: retrieve coordinates in dataframe form

require(deldir)

# Defining the coordinates 
x <- coor$X
y <- coor$Y

# Calculate the Delaunay triangulation, then the tiles.
z <- deldir(x,y)
w <- tile.list(z)


# Make a list of colours, and use 'em to plot:
ccc <- terrain.colors(18) 
plot(w,fillcol=ccc,close=TRUE)

library('dismo')

#Getting the coordinates
coor=as.data.frame(st_coordinates(bunnpris_sf))

#Converting x and y coordinates into a matrix
x <- coor$X
y <- coor$Y
points <- matrix(c(x,y), ncol=2)


#Implementing the triangulation and setting the boundaries as min and max of x and y coordinates of Bergen area
vor <- voronoi(points,ext = c(-40733.07,-10356.25,6714213,6749144) )


#plotting
spplot(vor, "id")

#Checking the class of the returned object 
class(vor)

[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"

5.2.Market Geo-segmentation

In this section, series of operations are implemented in order to calculate the total number of potential customer “population” in each tringulated area for each store.

Now, the triangulated areas can be plotted on Bergen map.

tmap_mode("view")
tm_shape(Bergen_popluation)+ tm_fill(col ="pop_tot" , alpha = .8 , n = 10) + tm_shape(bunnpris_sf) + tm_symbols(shape=1,size =.3,col = "red" )  + tm_shape(vor_sf) +tm_borders(col = "blue", lwd = 2, lty = "solid")+tm_text("id")

The following chuck of code join the triangulation polygon with Bergen_population data and aggregated the population on each triangulated area.

Bergen_popluation_tringulated = Bergen_popluation %>%
  st_join(vor_sf) %>%
  group_by(id) %>%
  summarize(total_population_for_each_store = sum(pop_tot, na.rm = TRUE))

Visualizing the distribution of potential market for the 18 locations.

On the left graph, it can be seen the total number of potential customers for each store labeled by its ID. The vertical blue line is the mean of potential customers across all the the stores. Assuming that the mean is the standard population number that each store should serve. We can see that there are some areas which are under-served such as; 14,15 and 18. Since these areas are not on the center of city, that gives an opportunity for establishing new stores. Areas on the city-center such as 1,2,3 and 4 can be considered as over-served. However, it has to be taken into consideration that there is inflow “mobility” of people to the center and potential customers are not only the people who are living on these areas, but also people who come to work or visit the center everyday and they are not represented on the data.